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Abstract 

The theory of monotone dynamical systems has been found very useful in the modeling of some 
gene, protein, and signaling networks. In monotone systems, every net feedback loop is positive. On 
the other hand, negative feedback loops are important features of many systems, since they are required 
for adaptation and precision. This paper shows that, provided that these negative loops act at a 
comparatively fast time scale, the main dynamical property of (strongly) monotone systems, convergence 
to steady states, is still valid. An application is worked out to a double-phosphorylation "futile cycle" 
motif which plays a central role in cukaryotic cell signaling. 

1 Introduction 

Monotone dynamical systems constitute a rich class of models, for which global and almost-global con- 
vergence properties can be established. They are particularly useful in biochemical applications and also 
appear in areas like coordination [28] and other problems in control [7j- One of the fundamental results 
in monotone systems theory is Hirsch's Generic Convergence Theorem [17 \ \18 \ \19 \ I37j. Informally stated, 
Hirsch's result says that almost every bounded solution of a strongly monotone system converges to the 
set of equilibria. There is a rich literature regarding the application of this powerful theorem, as well as 
of other results dealing with everywhere convergence when equilibria are unique ([9"1 122} 137]). to models of 
biochemical systems. See for instance [391 BO] f° r expositions and many references. 

Unfortunately, many models in biology are not monotone, at least with respect to any standard orthant 
order. This is because in monotone systems (with respect to orthant orders) every net feedback loop should 
be positive, but, on the other hand, in many systems negative feedback loops often appear as well, as they 
are required for adaptation and precision. However, intuitively, negative loops that act at a comparatively 
fast time scale should not affect the main characteristics of monotone behavior. The main purpose of this 
paper is to show that this is indeed the case, in the sense that singularly perturbed strongly monotone 
systems inherit generic convergence properties. A system that is not monotone may become monotone once 
that fast variables are replaced by their steady-state values. In order to prove a precise time-separation 
result, we employ tools from geometric singular perturbation theory. 

This point of view is of special interest in the context of biochemical systems; for example, Michaelis 
Menten kinetics are mathematically justified as singularly perturbed versions of mass action kinetics |11[ 
[29] . One particular example of great interest in view of current systems biology research is that of "futile 
cycle" motifs, as illustrated in Figure [TJ As discussed in |34j . futile cycles (with any number of intermediate 
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Figure 1: Dual futile cycle. A substrate So is ultimately converted into a product P, in an "activation" 
reaction triggered or facilitated by an enzyme E, and, conversely, P is transformed back (or "deactivated" ) 
into the original Sq, helped on by the action of a second enzyme F. 

steps, and also called substrate cycles, enzymatic cycles, or enzymatic interconversions) underlie signal- 
ing processes such as GTPase cycles [10], bacterial two-component systems and phosphorelays [31 [15] actin 
treadmilling [6]), and glucose mobilization |24j . as well as metabolic control [H] and cell division and apop- 
tosis [12] and cell-cycle checkpoint control [26]. A most important instance is that of Mitogen- Activated 
Protein Kinase ("MAPK") cascades, which regulate primary cellular activities such as proliferation, differ- 
entiation, and apoptosis [21 [201 US] i n eukaryotes from yeast to humans. MAPK cascades usually consist 
of three tiers of similar structures with multiple feedbacks [I] , [13] , [U] . Here we focus on one individual 
level of a MAPK cascade, which is a futile cycle as depicted in Figured! The precise mathematical model is 
described later. Numerical simulations of this model have suggested that the system may be monostable or 
bistable, see [27J. The later will give rise to switch-like behavior, which is ubiquitous in cellular pathways 
([14, 32, 35, 36J). In either case, the system under meaningful biological parameters shows convergence, 
not other dynamical properties such as periodic behavior or even chaotic behavior. Analytical studies done 
for the quasi-steady-state version of the model (slow dynamics), which is a monotone system, indicate that 
the reduced system is indeed monostable or bistable, see [HI]. Thus, it is of great interest to show that, 
at least in certain parameter ranges (as required by singular perturbation theory), the full system inherits 
convergence properties from the reduced system, and this is what we do as an application of our results. 

A feature of our approach, as in other control problems [I], [21], is the use of geometric invariant 
manifold theory |12l [23l 30J. There is a manifold M e , invariant for the full dynamics of a singularly 
perturbed system, which attracts all near-enough solutions. However, we need to exploit the full power of 
the theory, and especially the fibration structure and an asymptotic phase property. The system restricted 
to the invariant manifold M e is a regular perturbation of the slow (e=0) system. As remarked in Theorem 
1.2 in Hirsch's early paper [T7], a C 1 regular perturbation of a flow with eventually positive derivatives 
also has generic convergence properties. So, solutions in the manifold will generally be well-behaved, and 
asymptotic phase implies that solutions near M e track solutions in M £ , and hence also converge to equilibria 
if solutions on M e do. A key technical detail is to establish that the tracking solutions also start from the 
"good" set of initial conditions, for generic solutions of the large system. 

A preliminary version of these results was presented at the 2006 Conference on Decision and Control, 
and dealt with the special case of singularly perturbed systems of the form: 

x =f{x,y) 
ey =Ay + h(x) 

on a product domain, where A is a constant Hurwitz matrix and the reduced system x = f(x, —A~ 1 h(x)) 
is strongly monotone. However, for the application to the above futile cycle, there are two major problems 
with that formulation: first, the dynamics of the fast system have to be allowed to be nonlinear in y, and 
second, it is crucial to allow for an e-dependence on the right-hand side as well as to allow the domain to 
be a polytope depending on e. We provide a much more general formulation here. 

We note that no assumptions are imposed regarding global convergence of the reduced system, which 
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is essential because of the intended application to multi-stable systems. This seems to rule out the appli- 
cability of Lyapunov-theoretic and ISS tools [8l H3]. 

This paper is organized as follows. The main result is stated in Section [2j In Section [3l we review 
some basic definitions and theorems about monotone systems. The detailed proof of the main theorem 
can be found in Section [H and applications to the MAPK system and another set of ordinary differential 
equations are discussed in Section Finally, in Section [61 we summarize the key points of this paper. 

2 Statement of the Main Theorem 

In this paper, we focus on the dynamics of the following prototypical system in singularly perturbed form: 

-^ = f (x,y,e) (1) 
dy , . 

We will be interested in the dynamics of this system on a time- varying domain D £ . For < £ < 1, the 
variable x changes much slower than y. As long as e ^ 0, one may also change the time scale to r = t/e, 
and study the equivalent form: 

^ = efo(x,y,e) (2) 

dy , s 
~r = so(x,y,s). 

ar 

Within this general framework, we will make the following assumptions (some technical terms will be 
defined later), where the integer r > 1 and the positive number eo are fixed from now on: 

Al Let [/Cl n and V C M m be open and bounded. The functions 

fo:UxVx [0,e ] ^R n 

g :UxVx [0,e ]^IR m 

are both of class C£, where a function / is in C[ if it is in C r and its derivatives up to order r as 
well as / itself are bounded. 

A 2 There is a function 

m : U -> V 

of class C£, such that go(x, mo(x), 0) = for all x in U. 

It is often helpful to consider z = y — mo(x), and the fast system ([2]) in the new coordinates becomes: 

^-=ef 1 (x,z,£) (3) 
dz 

— = 31(2, z,e), 
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where 



fi(x,z,e) = f (x,z + m (x),e), 

gi(x,z,e) = g (x,z + m (x),e) - e[D a .m (x)]/i(x, z, e). 



When e = 0, the system 



(|3j) degenerates to 



— = g 1 (x,z,0), x(r)=x eU, 



(4) 



seen as equations on {z \ z + mo(xo) G V}. 

A3 The steady state z = of ([1]) is globally asymptotically stable on {z | z + mo(xo) G for all xo G £7. 

A4 All eigenvalues of the matrix D y go(x,mo(x),0) have negative real parts for every x G U, i.e. the 
matrix D y go(x, mo(x), 0) is Hurwitz on [/. 

A5 There exists a family of convex compact sets D e C U xV, which depend continuously on e G [0, £o]; 
such that (H|) is positively invariant on D e for e G (0,£o]. 

A6 The flow ifi® of the limiting system (set e = in ([TJ): 



onto the x-axis. 
A7 The set of equilibria of ([1]) on D e is totally disconnected. 

Remark 1 Assumption A3 implies that y = mo(x) is a unique solution of go(x, y, 0) = on U. 
Continuity in A5 is understood with respect to the Hausdorff metric. 

In mass-action chemical kinetics, the vector fields are polynomials. So, Al follows naturally. 
Our main theorem is: 

Theorem 1 Under assumptions Al to A7, there exists a positive constant e* < eq such that for each e G 
(0, £*), the forward trajectory of §\§ starting from almost every point in D e converges to some equilibrium. 

3 Monotone Systems of Ordinary Differential Equations 

In this section, we review several useful definitions and theorems regarding monotone systems. As we 
wish to provide results valid for arbitrary orders, not merely orthants, and some of these results, though 
well-known, are not readily available in a form needed for reference, we provide some technical proofs. 

Definition 1 A nonempty, closed set C C is a cone if 



— = / (x,m (x),0) 



(5) 



has eventually positive derivatives on Kq, where K$ is the projection of 
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1. C + C c C, 

2. R+C C C, 

3. Cf)(-C) = {0}. 

We always assume C ^ {0}. Associated to a cone C is a partial order on R . For any x,y G R , we 
define: 

x > y <±> x — y £ C 

x > y ^ x — y £ C, x y. 

When IntC is not empty, we can define 

i>i/<S>x-i/G IntC. 



Definition 2 T/ie duaZ cone o/C is defined as 

C* = {A G (R^)* | A(C) > 0}. 

An immediate consequence is 

x eC^X{x) > 0,VA g C* 
x G IntC ^ X(x) > 0, VA G C* \ {0}. 

With this partial ordering on R-^, we analyze certain features of the dynamics of an ordinary differential 
equation: 

where F : R N — > R^ is a C 1 vector field. We are interested in a special class of equations which preserve 
the ordering along the trajectories. For simplicity, the solutions of ([6]) are assumed to exist for alH > in 
the sets considered below. 

Definition 3 The flow <pt of ([6}) is said to have (eventually) positive derivatives on a set V C R^, if 
[D z <f> t (z)]x G IntC for all x G C \ {0}, z G V, and t>0 (t > t for some t > 0). 

It is worth noticing that [D z (j)t(z)]x G IntC is equivalent to \{\D z (f>t{z)\x) > for all A G C* with norm 
one. We will use this fact in the proof of the next lemma, which deals with "regular" perturbations in the 
dynamics. The proof is in the same spirit as in Theorem 1.2 of [18] . but generalized to the arbitrary cone 
C. 

Lemma 1 Assume V C R^ is a compact set in which the flow <fit of (0) has eventually positive derivatives. 
Then there exists 5 > with the following property. Let ipt denote the flow of a C 1 vector field G such that 
the C 1 norm of F{z) — G(z) is less than 5 for all z in V. Then there exists t* > such that if ip s (z) G V 
for all s G [0, t] where t > i*, then [D z ip t (z)]x G IntC for all z G V and x G C \ {0}. 
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Proof. Pick t > so that X([D z cj) t (z)]x) > for all t > t , z G V, A G C*, x G C with |A| = 1, \x\ = 1. 
Then there exists 5 > with the property that when the C 1 norm of -F(^) — G(z) is less than (5, we have 
\{[D z il) t {z))x) > for t < t < 2* . 

When t > 2io, we write t = r + fcto, where to < r < 2to and G N. If ip s {z) G F for all s G [0, t], we 
can define := ipjt (z) for j = 0, . . . , fc. For any x G C \ {0}, using the chain rule, we have: 

[D z ip t (z)]x = [D z ip r (zk)][D z ipt (zk-l)] ■ ■ ■ [D z ip tQ (zo)]x. 
By induction, it is easy to see that [D z ipt(zj\x G IntC. I 

Corollary 1 If V is positively invariant under the flow ipt, then ipt has eventually positive derivatives in 
V. 

Proof. If V is positively invariant under the flow ipt, then for any z G V the condition i/j s (z) G V for 
s G [0, t] is satisfied for all t > 0. By the previous lemma, tpt has eventually positive derivatives in V. I 

Definition 4 The system (|6|) or the flow <pt of (|6|) is called monotone (resp. strongly monotone) in a set 
W C R N , if for allt>0 and z u z 2 G W , 

Zl > Z 2 =>(pt( z l) > 4>t{z2) 

(resp. <f>t(zi) » 4>t( z 2) when z\ ^ z 2 ). 
It is eventually (strongly) monotone if there exists to > such that (ftt is (strongly) monotone for all t > to- 

Definition 5 An set W C M. is called p-convex, if W contains the entire line segment joining x and y 
whenever x < y, x, y G W . 

The next two propositions discuss the relations between the two definitions, (eventually) positive derivatives 
and (eventually) strongly monotone. 

Proposition 1 Let W C be p-convex. If the flow 4>t has (eventually) positive derivatives in W , then 
it is (eventually) strongly monotone in W . 

Proof. For any z\ > z 2 G W,X G C* \{0} and t > (t > to for some to > 0), we have that X((pt(zi) — 4>t(z 2 )) 
equals 

/ X([D z 4> t ( SZl + (1 - s)z 2 )){z 1 - z 2 ))ds > 0. 
Jo 

Therefore, (ftt is (eventually) strongly monotone in W. I 

Proposition 2 Suppose 4>t is (eventually) strongly monotone on an open set U C M. N . Then 4>t has 
(eventually) positive derivatives in U. 

Proof. Fix t > such that (f>t as a function from U to M N is strongly monotone (i.e. 4>t{z\) ^> <pt{z 2 ), 
whenever z\ > z 2 , z\ / z 2 in U). 

For any A G C* \ {0}, x G C \ {0}, and z G U, 

X(\D z (j) t (z)}x) = lim h^XiMz + hx) - (f> t (z))- 

h— >0 

If h > is small enough such that z+hx G U, then we have X((j)t{z-\-hx) — <ftt(z)) > 0. Thus X([D z (j)t(z)]x) > 
0, and [D z <p t (z)]x G IntC. I 
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Lemma 2 Suppose that the flow of (|6|) has compact closure and eventually positive derivatives in a 
p-convex set W C R . If the set of equilibria is totally disconnected (e.g. countable), then the forward 
trajectory starting from almost every point in W converges to an equilibrium. 

Proof. By Proposition [IJ <t>t is eventually strongly monotone. The result easily follows from Hirsch's 
Generic Convergence Theorem (|19j. |37j). I 



4 Details of the Proof 

Our approach to solve the varying domain problem is motivated by Nipp |30j . The idea is to extend the 
vector field from U x V x [0, e] to R n x M m x [0, £o], then apply geometric singular perturbation theorems 
([33J) on M n x M. m x [0, £q], and finally restrict the flows to D £ for the generic convergence result. 



4.1 Extensions of the vector field 

For a given compact set K C W 1 (Kq C K C U), the following procedure is adopted from [3D] to extend 
a CT function with respect to the x coordinate from U to W 1 , such that the extended function is C£ and 
agrees with the old one on K. This is a routine "smooth patching" argument. 

Let U\ be an open subset of U with C r boundary and such that K C U\ C U. For 0o > sufficiently 
small, define 

Uf° := {x G U\ | @(x) > ©o}, where Q(x) := min \x — u\, 

u£dUi 

such that K is contained in U®° . Consider the scalar C°° function p: 

a < 
p(a) := ^ exp(l — exp(a — l)/a) < a < 1 

1 a > 1. 



Define 



and 



x e R n \ Ui 
Q(x) := { Q(x) xeUt\ U®° 



Q x G C/® 



For any q^C r b {U), let 



f/(.c) : - ^ ^ | ^ ^ and : = <d(x)q(x) 



Then g(x) G C fe r (M n ) and q(x) = q{x) on K. 
We fix some c?o > such that 



D do := {z G R m | |z| < d } C P| {z I z + m (x) G F}. 



xeK 
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Then we extend the functions fi and mo to f\ and mo respectively with respect to x in the above way. To 
extend gi, let us first rewrite the differential equation for z as: 

dz 

— = [B(x) + C(x,z)]z + eH(x,z,e) - e[D x m (x)]fi(x, z, e), 

where 

B(x) = D y go(x,mo(x),0) and C(x,0) = 0. 

Following the above procedures, we extend the functions C and H to C and H , but the extension of B is 
defined as 

B(x) := G(x)B(x) - - G(x))/ n , 

where /U is the positive constant such that the real parts of all eigenvalues of B{x) is less than —\i for every 
x € K. According to the definition of B(x), all eigenvalues of B(x) will have negative real parts less than 
— ji for every x € W 1 . The extension gi, defined as: 

[B(x) +C(x,z)]z + eH(x,z,e) - e[D x m (x)]fi(x, z, e), 

is then C£ -1 (M ra x D do x [0,£o]) and agrees with g\ on K x D do x [0,£o]. 

To extend functions f\ and g\ in the z direction from Dd Q to M m , we use the same extension technique 
but with respect to z. Let us denote the extensions of f\,C,H and the function z = z by fx, C, H and z 
respectively, then define g± as: 

[B(x) + C(x,z)]z(z) +eH(x,z,e) - e[D x fh (x)]fi(x, z, e), 

which is now C^ _1 (M n x R m x [0, £o]) and agrees with g\ on K x Z)rf 1 x [0, Eq] for some d\ slightly less than 
do. Notice that z = is a solution of 5i(x,z,0) = 0, which guarantees that for the extended system in 
(x,y) coordinates (y = z + mo(x)): 

^ = ef(x,y,e) (7) 
= 9[x,y,e), 

y = fho(x) is the solution of g(x,y,0) = 0. To summarize, (J7|) satisfies 
El The functions 

xM m x [0,e ]), 
gE^rxR-x^o]), 
m G C fe r (M n ), 5 (x,m (a;),0) =0, VxG R n . 

E2 All eigenvalues of the matrix D y g(x, mo(x),0) have negative real parts less than —\x for every x £ R n . 

E3 The function mo coincides with mo on K, and the functions / and g coincide with /o and go respectively 
on 

^di ■= {(x,y) \ x G K, \y - m (x)\ < di}. 

Conditions El and E2 are the assumptions for geometric singular perturbation theorems, and condition 
E3 ensures that on the flow of ([2]) coincides with the flow of ©. If we apply geometric singular 
perturbation theorems to ([7]) on M. n x M. m x [0, £o] > the exact same results are true for ([2]) on fi^ . For the 
rest of the paper, we identify the flow of ([7]) and the flow of (|2]) on £1^1 without further mentioning this 
fact. (Later, in Lemmas UH7J when globalizing the results, we consider again the original system.) 
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4.2 Geometric singular perturbation theory 

The theory of geometric singular perturbation can be traced back to the work of Fenichel [12], which first 
revealed the geometric aspects of singular perturbation problems. Later on, the works by Knobloch and 
Aulbach [25j, Nipp [30J, and Sakamoto [33J also presented results similar to [12]. By now, the theory 
is fairly standard, and there have been enormous applications to traveling waves of partial differential 
equations, see |23] and the references there. For control theoretic applications, see [Tj 121]. 

To apply geometric singular perturbation theorems to the vector field on M. n x R m x [0, £q], we use the 
theorems stated in [33]. The following lemma is a restatement of the theorems in [33], and we refer to [33] 
for the proof. 

Lemma 3 Under conditions El and E2, there exists a positive E\ < Eq such that for every e G (0, £i]: 

1. There is a C^ -1 function 

m : R n x [0,£i] -> M. m 

such that the set M £ defined by 

M £ := {(x,m(x,e)) \ x G M. n } 
is invariant under the flow generated by Moreover, 

sup \m(x, e) — fh (x)\ = 0(e), as e — > 0. 

In particular, we have m(x,0) = fho(x) for all x G M. n . 

2. The set consisting of all the points (xo,yo) such that 

sup \y(T;x ,y ) - m(x(r; x , y ), e) \e%~ < oo, 

T>0 

where (x(T;xo,yo),y(T;xo,yo)) is the solution of ([?]) passing through (xo,yo) at t = 0, is a C r_1 - 
immersed submanifold in W l x R m of dimension n + m, denoted by W s (M e ). 

3. There is a positive constant 5q such that if 

sup |y(r; x ,y ) - m(x(r; x , yo),e) \ < 5 , 

T>0 

then (x ,y ) G W S (M £ ). 

4- The manifold W S (M £ ) is a disjoint union of C r_1 -immersed manifolds W^(£) of dimension m: 

W S {M £ ) = (J W £ s (0- 

For each t; G W 1 , let H £ (^)(t) be the solution for r > of 

— =ef{x,m(x,e),e), x(0) = f G M n . 
Then, the manifold W £ (£) is the set 

flT (AT 

{(x ,yo) I sup \x(r)\e~ < oo,sup \y(r)\e~ < oo}, 

T>0 T>0 
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where 

x( T ) = x(T;x ,y ) - fig (£)(>), 

y( T ) = y(r;x ,yo) - m(H £ (£)(r), e) . 

5. The fibers are "positively invariant" in the sense that W*(H e (£)(r)) is the set 

{(x(T;x ,yo),y(T;x ,yo)) I (x ,yo) e W £ (£)} 
for each r > 0, see Figured 

6. The fibers restricted to the 5q neighborhood of M £ , denoted by , can be parametrized as follows. 
There are two C£ _1 functions 

P eA) : R n x D 5o -> R n 
Q eM : M m x D So - R m , 

and a map 

T £A : M n x D 5o -.K"xr 

mapping (£,77) to (x,y), where 

x = i + P £ ,s (£,77), y = m(x, e) + Q £)5o (£, 77) 

w/ A (0 = r eA (e,^). 

■ 

Remark 2 TTie <5o Jtj property 3 can be chosen uniformly for e € (0, £0] . Without loss of generality, we 
assume that 5q < d\. 

Notice that property 4 insures that for each (xo,ya) € W S {M £ ), there exists a £ such that 

\x(r;x ,yo) ~ H e (£)(t)\ -> 0, 
\y{T;x ,y ) - m(H £ (£)(T),e)\ -> 0. 
as t — > 0. Taws is q/ton referred as the "asymptotic phase" property, see Figure® 

4.3 Further analysis of the dynamics 

The first property of Lemma [3] concludes the existence of an invariant manifold M £ . There are two reasons 
to introduce M £ . First, on M £ the ^-equation is decoupled from the y-equation: 

^ = f(x,m{x,e),e) (8) 
y(t) = m(x(t),e). 

This reduction allows us to analyze a lower dimensional system, whose dynamics may have been well 
studied. Second, when e approaches zero, the limit of (jSJ) is ([5]) on Kq. If © has some desirable property, 
it is natural to expect that this property is inherited by ([8]). An example of this principle is provided by 
the following Lemma: 
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Figure 2: An illustration of the "positive invariant" and "asymptotic phase" properties. Let po be a point 
on the fiber W £ {qo) (vertical curve). Suppose the solution of (|7|) starting from qo € M £ evolves to q\ £ M £ 
at time n, then the solution of starting from pq will evolve to p\ E W £ (q\) at time t\. At time T2, they 
evolve to (?2j2?2 respectively. These two solutions are always on the same fiber. If we know that the one 
starting from qo converges to a equilibrium, then the one starting from pq also converges to a equilibrium. 

Lemma 4 There exists a positive constant e<i < £\, such that for each e G (0,62), the flow ?/>f of ([8]) has 
eventually positive derivatives on K £ , which is the projection of M £ f^\D £ to the x-axis. 

Proof. Assumption A6 states that the flow ifi® of the limiting system ([5|) has eventually positive derivatives 
on Kq. By the continuity of m{x, e) and D £ at e=0, we can pick £2 small enough such that the flow ip® has 
eventually positive derivatives on K £ for all e € (0,62)- Applying Corollary [Q we conclude that the flow 
V'f °f ® has eventually positive derivatives on K £ provided K £ is positively invariant under dH), which 
follows easily from the facts that ([7]) is positively invariant on D £ and M £ is an invariant manifold. I 

The next Lemma asserts that the generic convergence property is preserved for ([8]), see Figure [3j 



Figure 3: This is a sketch of the manifolds Mq (surface bounded by dashed curves), M £ (surface bounded 
by dotted curves), and D £ (the cube). It highlights two major characters of M £ . First, M £ is close to Mq. 
Second, the trajectories on M £ converge to equilibria if those on Mq do. 

Lemma 5 For each e G (0, 82), there exists a set C £ C K £ such that the forward trajectory of ([8]) starting 
from any point of C £ converges to some equilibrium, and the Lebesgue measure of K £ \ C £ is zero. 
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Proof. Apply Lemma [2] and Lemma E] under assumptions A5 and A7. 



By now, we have discussed flows restricted to the invariant manifold M £ . Next, we will explore the 
conditions for a point to be on W S (M £ ), the stable manifold of M £ . Property 3 of Lemma [3] provides a 
sufficient condition, namely, any point (xo,yo) such that 

sup|y(r;x ,y ) - m(a?(r; x , yo), e) \ < S (9) 

r>0 

is on W s (M e ). In fact, if we know that the difference between yo an d m(xQ,e) is sufficiently small, then 
the above condition is always satisfied. More precisely, we have: 



Lemma 6 There exists £3 > 0, 5q > d > 0, such that for each e E (0,63), if the initial condition satisfies 
I yo —m(xo,e)\ < d, then ([9]) holds, i.e. (xo,yo) G W S (M £ ). 



Proof. Follows from the proof of Claim 1 in [30]. I 

Before we get further into the technical details, let us give an outline of the proof of the main theorem. 
The proof can be decomposed into three steps. First, we show that almost every trajectory on D £ f] M £ 
converges to some equilibrium. This is precisely Lemma 03 Second, we show that almost every trajectory 
starting from W S {M £ ) converges to some equilibrium. This follows from Lemma [5] and the "asymptotic 
phase" property in Lemma [31 but we still need to show that the set of non-convergent initial conditions is 
of measure zero. The last step is to show that all trajectories in D £ will eventually stay in W S (M £ ), which 
is our next lemma: 



Lemma 7 There exist positive tq and £4 < £3, such that (x(to), 2/(to)) € W S (M £ ) for all e € (0, £4), where 
(x(r),y(r)) is the solution to ([2]) with the initial condition (xo,yo) G D E . 

Proof. It is convenient to consider the problem in (x,z) coordinates. Let (x(t),z(t)) be the solution to 
([3]) with initial condition (xq, zq), where zq = yo — m(xo, 0). We first show that there exists a tq such that 
|z(t )| <d/2. 

Expanding gi(x, z, e) at the point (xq,z, 0), the equation of z becomes 

dz dgi . 

j- = 9i{xo, z, 0) + -Q—KS, z, 0){x - x ) + eR{x, z, e) 
for some ^(r) between xq and x(t) (where ^(r) can be picked continuously in r). Let us write 

z{t) = z°(t) + w(t), 

where z°(t) is the solution to (|1J) with initial the condition z°(0) = zq, and w(t) satisfies 

— = gi(x ,z,0) - gi(x ,z°,0) + —(£,z,Q)(x - x ) + eR(x, z, e) (10) 

= ^( x o, C, 0)w + z, 0) J fi(x(s), z(s),e) ds + eR(x, z, e), 

with the initial condition w(0) = and some C( T ) between and z(t) (where £(r) can be picked 

continuously in r). 

By assumption A3, there exist a positive tq such that |-z°(r)| < d/A for all r > To. Notice that we are 
working on the compact set D £ , so tq can be chosen uniformly for all initial conditions in D £ . 
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We write the solution of (jlOp as: 

w(t)=J ^-(x ,(,0)wds + e J (^-(£,z,0) J f x (x, z,s) ds' + R(x, z,e)) ds. 
Since the functions f±, R and the derivatives of g\ are bounded on D £ , we have: 

\w(t)\ < J L\w\ds + e J J M 2 ds' + M^j ds, 

for some positive constants L, Mj, z = 1,2,3. The notation means the Euclidean norm of w G K m . 
Moreover, if we define 

a(r) = jf ^Mi M 2 ds' + M 3 ^j ds, 

then 

|iw(r)| < / ds + ea(ro), 

for all r G [0, ro] as a is increasing in r. Applying Gronwall's inequality ([38]), we have: 

\w(t)\ < ea(ro)e Lr , 

which holds in particular at r = ro. Finally, we choose £4 small enough such that ea(To)e LT ° < cZ/4 and 
\m(x,e) — m(x, 0)| < d/2 for all e G (0,64). Then we have: 

\y(ro) - m(x(T ),e)\ < |y(r ) - m(x(r ),0)[ + |m(x(r ),e) - m(x(r ),0)| 

< \z(t )\ + d/2 

< d/2 + d/2 = d. 

That is, (x(ro),y(r )) G I<F s (M e ) by Lemma I 
By now, we have completed all these three steps, and are ready to prove Theorem [IJ 

4.4 Proof of Theorem [1] 

Proof. Let e* = min{e2,£4}- For e G (0, e*), it is equivalent to prove the result for the fast system ([2]). 
Pick an arbitrary point (xo,yo) m an< ^ t nere are three cases: 

1. yo = m(xo,e), that is, (xo,yo) G M e P|D e . By Lemma the forward trajectory converges to an 
equilibrium except for a set of measure zero. 

2. < \yo — m(xo,s)\ < d. By Lemma El we know that (xo,yo) is in W S (M £ ). Then, property 4 of 
Lemma [3] guarantees that the point (xo,yo) is on some fiber W* where £ G K £ . If £ G C £ , that is, 
the forward trajectory of £ converges to some equilibrium, then by the "asymptotic phase" property 
of Lemma [3l the forward trajectory of (xo>Z/o) a ^ so converges to an equilibrium. To deal with the 
case when £ is not in C £ , it is enough to show that the set 

Be,d= U W ld(0 
t<BK e \C e 
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has measure zero in W m+n . Define 

S £ ,d = (K £ \ C £ ) x D d . 

By Lemma O K £ \ C £ has measure zero in M. n , thus S £t d has measure zero in R™ x M. m . On the other 
hand, property 6 in Lemma [3] implies B £ d = T Ei d(S Et d)- Since Lipschitz maps send measure zero sets 
to measure zero sets, B £ d is of measure zero. 

3. |yo — r n(xQ,e)\ > d. By LemmaEl the point (x(to), v(to)) is in W S (M £ ) and we are back to case 2. 
The proof is completed if the set 4> e _ T(j (B £) d) has measure zero, where (f> e T is the flow of ([2]). This is 
true because (j) e T is a diffeomorphism for any finite r. 



5 Applications 

5.1 An application to the dual futile cycle 

We assume that the reactions in Figure [T] follow the usual enzymatic mechanism (|13|): 

k\ k 3 
S + E^d H S x + E^C 2 h S 2 + E 

h\ h 3 

S 2 + F^C 3 H Si + F^C 4 H So + F. 
h-i h- 3 

There are three conservation relations: 

Stot = [So] + [Si] + [S 2 ] + [Ci] + [C 2 ] + [C 4 ] + [C 3 ], 
E tot = [E] + [C 1 ] + [C 2 ], 
Ftot = [F] + [C 4 ] + [C 3 ], 

where brackets indicate concentrations. Based on mass action kinetics, we have the following set of ordinary 
differential equations: 



d[S Q ] 
dr 
d[S 2 ] 

dr 
d[Ci] 

dr 

d[C 2 ] 

dr 
d[C 4 ] 

dr 
d[C 3 ] 

dT 



h 4 [C4-ki[So][E]+k_i[Ci] 
h[C 2 ] - hi[S 2 ][F] + h_i[C 3 ] 

ki[S ][E} - (k^i + k 2 )[Ci] (11) 
k 3 [Si}[E] - (k_ 3 + k 4 )[C 2 ] 
h 3 [Si}lF] - (h_ 3 + h 4 )[C4} 
h 1 [S 2 ][F]-(h- 1 + h 2 )[C 3 ]. 
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After rescaling the concentrations and time, (fTT|) becomes: 
—— = -k\StatXx{l -yi- y 2 ) + fe— 12/1 + ^2/3 

-7T = -hiS to tcx 2 (l - 2/3 - 2/4) + fo-iq/4 + ^42/2 
at 

e-^- = kiS to txi(l - 2/1 - 2/2) - (fe-i + £2)2/1 (12) 
£ -JT = faStoti 1 -xi-x 2 - eyi - £2/2 - ecy 3 - ecy 4 ) x (1 - 2/1 - 2/2) - (k-3 + £4)2/2 



where 



= hsStatO- -xi-x 2 - eyi - ey 2 - ecy 3 - ecy 4 ) x (1 - y 3 - 2/4) - (h- 3 + h 4 )y 3 

at 

= hiStatX2(l - 2/3 - 2/4) - (h-i + h 2 )y 4 , 



[So] [S2} [Ci] [C 2 

xi = — , x 2 = 75 -) 2/1 = V 2 . 2/2 - 



Sfot "Siot ^tot -S'tot 

[C4] [C 3 ] £^ ot F iot 

2/3 = -^—, 2/4 = ^—, e = — , c=— , i = re. 

rtot r t ot otot &tot 

These equations are in the form of (pQ). The conservation laws suggest taking £0 = 1/(1 + c ) an d 

D £ ={(xi,x 2 ,2/i,2/2,2/3,2/4) |0 < 2/1 + 2/2 < 1,0 < 2/3 + 2/4 < 1, 
xi > 0, x 2 > 0, < xi + x 2 + £(2/1 + 2/2 + C2/3 + C2/4) < 1}- 

For £ £ (0, £0], taking the inner product of the normal of dD £ and the vector field, it is easy to check that 
(|12[) is positively invariant on D e , so A5 holds. We want to emphasize that in this example the domain 
D e is a convex polytope varying with e. 

It can be proved that on D e system (12) has at most a finite number of steady states, and thus A7 holds. 
This is a consequence of a more general result, proved using some of the ideas given in [16], concerning the 
number of steady states of more general systems of phosphorylation/dephosphorylation reactions, see 

Solving 27 (x,2/,0) = 0, we get 



yi 



y% 



3/3 



2/4 



K m i _|_ K m \(\-x\-xi) 

Stot K m 2 



Kml (1— gl— gg) 



K„ 



K m l _j Kml(l-X 1 -X 2 ) , ' 

Stot KmS 

^m4 

Km3 I ^3(1-11-12) 1 ^ 
Stot #m 4 
^2 

St„t X m4 "^ X2 



where K m i, K m2 , K m3 and i^ m 4 are the Michaelis-Menten constants defined as 

r , fe_i + A: 2 r . A:_ 3 + fc 4 /i-i + n 2 „ /i-3 + hi 

J^ml — ; , A m2 — : , A m 3 — , A m 4 — . 

fei fa fi1 Al3 
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Now, we need to find a proper set U C M? satisfying assumptions A1-A4. Suppose that U has the 
form 

U = {(xi, x 2 ) | x\ > -a, x 2 > -a, x\ + x 2 < 1 + a}, 

for some positive a, and V is any bounded open set such that D e is contained in U x V, then Al follows 
naturally. Moreover, if 



a < <jq : = min 



K m iK m2 K m ^K m ^ 



Stot{K m i + K m2 ) ' St t(K m 3 + K m 4) J 
A2 also holds. To check A4, let us look at the matrix: 



B{x) := D y g (x, m (x), 0) 

where 



Si(x) 
B 2 (x) J ' 



B ( x \ = I -hStotxi - [k-i + k 2 ) -kiStotXi 
~ )] ' 1 -hS to t{l - xi - x 2 ) -kzS to t{l - %\ - x 2 ) - (k- 3 + k±) /' 



and 

-hzStotil - xi - x 2 ) - (h- 3 + hi) -hzStotil -x\ - x 2 ) 
-hxSt ot x 2 -h\StotX 2 - (h-i + h 2 



B 2 (x) 



If both matrices B\ and B 2 have negative traces and positive determinants, then A4 holds. 
Let us consider B\ first. The trace of B\ is 

-kxStotxi - (k-x + k 2 ) - k 3 S tot (l -x\- x 2 ) - (/c_ 3 + /c 4 ). 

It is negative provided that 

k-x + k 2 + k- 3 + ki 
Stat(kx + k 3 ) 

The determinant of Bx is 

kx(k- 3 + k 4 )S to tXx + A: 3 (A;_i + fe2)iStot(l - Si - x 2 ) + (k-x + k 2 )(k- 3 + k 4 ). 

It is positive if 

(k-x + k 2 )(k- 3 + h) 

a < 



Stot (kx(k-3 + k A ) + k 3 (k-x+k 2 ))' 
The condition for B 2 can be derived similarly. To summarize, if we take 

^-1 + ^2 + ^-3 + ^4 (k-x + k 2 )(k-3 + ki) 



a = mm < uq 



5*at(fci + k 3 ) ' 5 to i(A:i(A:_ 3 + fc 4 ) + k 3 (k-x + fc 2 )) 
h-x + h 2 + h-3 + h 4 (h-x + h 2 )(h-3 + hi) 



S t ot(hx + /13) ' S to t(hx(h-3 + /i 4 ) + ^s(/i-i + h 2 )) 

then the assumptions Al, A2 and A4 will hold. 

Notice that y in (|12p is linear in y when e = 0, so gx (defined as in ([3])) is linear in z, and hence the 
equation for z can be written as: 

dz 

— = B(x )z, xq € U, 
where the matrix B(xq) is Hurwitz for every xq E U. Therefore, A3 also holds. 
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To check A6, let us look at the reduced system (e = in (|12|)): 

dxi hixi ft 4C Km4 

dt K m , | K ml (l- Xl -x 2 ) + AT m3 ^(l-ax-xa) . " ^H^'^J U^j 

Stat ^ K m2 S iot ^ K m4 

^ = b^i + h ^ . = F2(X1 X2) 

Stot ^ K mi "^ X 2 Stot f Km2 tJ-i 

It is easy to see that F\ is strictly decreasing in X2, and F2 is strictly decreasing in x\ on 

-Ko = {(^1,^2) I %i > 0, 2:2 > 0,xi + x 2 < 1}. 
So, (|13p is strongly monotone on some open set W containing Kq with respect to the cone 

{(xi,x 2 ) I xi < 0,x 2 > 0}. 

Applying Lemma [21 the flow of (|13p has eventually positive derivatives on Kq (c W), and A6 is valid. 



So the system formulated in the form of (112j) satisfies all assumptions Al to A7. Applying Theorem 
[H we have: 



Theorem 2 There exist a positive e* < £q such that for each e € (0,e*) ; the forward trajectory of (|12p 
starting from almost every point in D £ converges to some equilibrium. 

In fact, since the reduced system is of dimension two, we know that every trajectory in D e , instead of 
almost every trajectory in D e , converges to some equilibrium (|19j). 

It is worth pointing out that the conclusion we obtained from the above theorem is only valid for 
small enough e; that is, the concentration of the enzyme should be much smaller than the concentration 
of the substrate. Unfortunately, this is not always true in biological systems, especially when feedbacks 
are present. However, if the sum of the Michaelis-Menten constants and the total concentration of the 
substrate are much larger than the concentration of enzyme, a different scaling: 

[So] [S2] , E t ot , / 

XI = -J", X 2 = -j", £ = — , t = TE, 

where A = Stat + Kmi + Km2 + ^m3 + will allow us to obtain the same convergence result. 
5.2 Another example 

The following example demonstrates the importance of the smallness of e. Consider an m + 1 dimensional 
system: 

dx 

— = 7 (yi,...,y m ) (14) 

e-J 1 = - diyi- <Xi(x), di>0, i = l,...,m. 

under the following assumptions: 

1. There exists an integer r > 1 such that the derivatives of 7,/3, and a« are of class C£ for sufficiently 
large bounded sets. 
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2. The function (3{x) is odd, and it approaches infinity as x approaches infinity. 

3. The function cti(x) (i = 1, . . . , m) is bounded by positive constant Mj for all i£R. 

4. The number of roots to the equation 

7(ai(x), . . .,a m (x)) = (3(x) 

is countable. 

We are going to show that on any large enough region, and provided that e is sufficiently small, almost 
every trajectory converges to an equilibrium. To emphasize the need for small e, we also show that when 
e > 1, limit cycles may appear. 

Assumption 4 implies A7, and because of the form of (|14p . A3 and A4 follow naturally. A6 also holds, 
as every one dimensional system is strongly monotone. For A5, we take 

D e = {(x,y) \ \x\ < a, \yi\ <bi,i = l,... ,m}, 

where 6, is an arbitrary positive number greater than ^ and a can be any positive number such that 

/3(a) > N b := max 7(7/1, ... ,y m )- 

M<bi 

Picking such hi and a assures 

dx dyi 

< °> Vi-JT < °' 
dt dt 

i.e. the vector field points transversely inside on the boundary of D £ . Let U and V be some bounded open 
sets such that D £ C U x V, and assumption 1 holds on U and V. Then Al and A2 follow naturally. By 
our main theorem, for sufficiently small e, the forward trajectory of (|14p starting from almost every point 
in D e converges to some equilibrium. 

On the other hand, convergence does not hold for large e. Let 

(5(x) = — - x, ai(x) = 2tanhx, m = 1, 7(7/1) = y 1} d\ = 1. 

It is easy to verify that (0, 0) is the only equilibrium, and the Jacobian matrix at (0, 0) is 

1 1 

-2/e -1/e 

When £ > 1, the trace of the above matrix is 1 — 1/e > 0, its determinant is 1/e > 0, so the (only) 
equilibrium in D is repelling. By the Poincare-Bendixson Theorem, there exists a limit cycle in D. 



6 Conclusions 

Singular perturbation techniques are routinely used in the analysis of biological systems. The geometric 
approach is a powerful tool for global analysis, since it permits one to study the behavior for finite e 
on a manifold in which the dynamics is "close" to the slow dynamics. Moreover, and most relevant to 
us, a suitable fibration structure allows the "tracking" of trajectories and hence the lifting to the full 
system of the exceptional set of non-convergent trajectories, if the slow system satisfies the conditions of 
Hirsch's Theorem. Using the geometric approach, we were able to provide a global convergence theorem 
for singularly perturbed strongly monotone systems, in a form that makes it applicable to the study of 
double futile cycles. 
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